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Abstract: Coupled arrays of Andronov-Hopf oscillators are investigated. These arrays can be 
diffusively or repulsively coupled, and can serve as central pattern generator models in animal 
locomotion and robotics. It is shown that repulsive coupling generates out-of-phasc oscillations, 
while diffusive coupling generates synchronous oscillations. Specifically, symmetric solutions 
and their corresponding amplitudes are derived, and contraction analysis is used to prove 
global stability and convergence of oscillations to either symmetric out-of-phase or synchronous 
states, depending on the coupling constant. Next, the two mechanisms are used jointly by 
coupling multiple arrays. The resulting dynamics is analyzed, in a model inspired by the CPG- 
motorneuron network that controls the heartbeat of a medicinal leech. 



1. I. INTRODUCTION 

Central pattern generators (CPGs) are often modeled as 
coupled nonlinear oscillators delivering phase-locked sig- 
nals. Some of their applications include animal locomo- 
tion, robotics (see Seo (2007), Ijspeert (2007), and Ijspeert 
(2008)), and other biological rhythmic behaviors such as 
e.g. the generation of a heartbeat, Buono (2004). In this 
paper, we explore a system of coupled Andronov-Hopf 
oscillators that can be used to model such phenomena. Re- 
pulsive coupling between the oscillators creates a traveling 
wave that simulates salamander gate (Ijspeert (2007)), and 
with additional coupling architecture, the heartbeat of a 
leech (Buono (2004)). 

Explicit rotational coupling between neighbouring oscilla- 
tors has been used in earlier work (sec Pham (2007), Seo 
(2007)) to generate a traveling wave or out-of-phase state. 
Here, we use inhibitory (repulsive) coupling to achieve a 
similar effect in perhaps a more physical way, allowing 
indeed the system itself to compute couplings achieving 
the out-of-phase behavior. While all our derivations aim 
at establishing global convergence results, we found some 
significant qualitative differences between these two types 
of couplings, such as the existence of two degenerate out- 
of-phase states in repulsive coupling, but not in rotational 
coupling. In addition, while rotational coupling preserves 
the amplitude of the coupled oscillators, repulsive coupling 
increases the amplitude of oscillation. In fact, this increase 
in amplitude is the key factor behind generating high 
frequency oscillation patterns in the heartbeat of a leech 
model, to be analyzed in the paper. 

The paper is organized as follows. In Part II, we ob- 
tain out-of-phase and synchronous solutions for nearest 
neighbor coupled Andronov-Hopf oscillators, and derive 
the condition for the existence of an out-of-phase state. 
In Part HI, either the out-of-phase or the synchronous 



solution is shown to be globally stable using techniques 
from contraction theory (for original derivation of contrac- 
tion theory see Lohmiller (1998)). The stability of the two 
different types of behavior is determined by the values of 
the coupling constant, with a positive coupling constant 
generating the out-of-phasc state and a negative constant 
resulting in a globally stable synchronous solution. In 
Part IV, the synchronous and the out-of-phase arrays are 
globally coupled to each other in a model inspired by the 
CPG in the heartbeat of a leech. Building on prior results, 
the bifurcation value for the onset of high frequency oscil- 
lations in the synchronous array is derived. 

2. II. STEADY-STATE DYNAMICS 

In this section we solve for the steady-state dynamics of 
the nearest neighbor coupled Andronov-Hopf oscillators, 
showing that both synchronous and out-phase solutions 
exist. In the next section, we consider the global stability 
properties of these solutions applying techniques from 
contraction theory. We first consider a ring of nearest 
neighbor coupled Andronov-Hopf oscillators, of the form: 

Xj = F(xj) -|-fc(xj -Xj+i -Xj_i) (1) 

where xj = {xj, yj} is a two dimensional vector describing 
the dynamics of the jth oscillator, with F(x) given by: 
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In complex form, Eqns. (1) and (2) can be expressed as: 

Zj = {a + iuj)zj - \zj fz.j + k (zj - - z^+i) (3) 

where z is a complex variable, given by: z = x ^iy. 

In the absence of coupling, the dynamics are that of a 
limit cycle, of amplitude, |xj| = ^Ja. and frequency w. 




Fig. 1. Repulsively coupled array, fc > 0, for even number 
of oscillators, N — 4. Shows that the neighboring 
oscillators are maximally out-of-phase (by tt) 

Due to symmetry considerations, the synchronous state 
is one possible solution to the above equation, re sulting 
in the amplitude of oscillation given by: |xj| = \/a — k. 
This solution will be shown in the next section to be 
globally stable for diffusive type of coupling, given by: 
A; < 0. For the "repulsive couphng", given by fc > 0, 
the system tends to an out-of-phase state, whereby the 
neighboring oscillators are maximally out of phase with 
each other. Landsman (2006). For oscillators coupled in 
a ring, this results in two different types of dynamics, 
depending on whether the number of oscillators, TV is 
even or odd. When N is even, the array oscillates with 
a difference of tt between nearest neighbors, splitting into 
two equal synchronous groups, that are 180 degrees out of 
phase with each other (see Figure 1, plotted for = 4). 
The phase difference between nearest neighbors is thereby 
given by: 

A</.— = TT (4) 

In the case of oscillators in a line, coupled with fc > 0, the 
above phase difference is the only globally stable solution. 
The situation becomes more complicated in a ring coupled 
model when N is odd. In this case the phase difference 
of TT between nearest neighbors is not a symmetric or a 
stable solution. For example, imagine a ring of 3 oscillators, 
where the 2nd oscillator is out of phase with the 1st and 
the 3rd by tt. Then the 1st and the 3rd oscillator will 
actually be in-phase, which is an unstable state. Requiring 
the neighboring oscillators to be maximally out of phase, 
while preserving the symmetry of the system leads to two 
possible degenerate solutions given by: 

A0°f+i-7r±7r/iV (5) 

The smallest phase difference between the two oscillators 
is given by the next to nearest neighbor phase difference: 
^(f>°f+2 ± 2Tr/N. The vectors xj therefore fall on a circle 
where they are spaced with an equal phase difference of 
2t:/N, forming a symmetric out-of-phase solution. 

Figure 2 illustrates the stable out-of-phase dynamics for 
N — 5. Each oscillator in the Figure is shifted in phase 



Fig. 2. Repulsively coupled array for odd number of 
oscillators, = 5. The neighboring oscillators are 
out of phase by n + 7t/5 (for example, compare the 
first and the second oscillator). 

from its neighbor by ir + ir/b. Comparing Eq. (5) to 
Eq. (4), we can see that as iV — >■ c», the steady-state 
dynamics of the A^-odd array approach that of iV-even. 
This should be expected, since for large N, adding one 
more oscillator to the ring (and thereby changing N from 
odd to even or vice versa) will not significantly change 
the energy function, which is minimized when Eqs. (5) 
and (4) for iV-odd and A^-even, respectively, are satisfied. 
To solve for the amplitude of the out-of-phase solution, 
which is stable for repulsive coupling, we use Eqn. (5), 
writing {zj-i = e±^('^+'^/^)zj,Zj+i = e^^^'^+'^/^^Zj} and 
substituting for {zj„i,Zj+i} into Eqn. (3). After grouping 
the linear terms, we can now solve for the amplitude of 
oscillation, getting: 

\z^\N-odd ^ + ^ (1 + 2cos{Tr/N))]^/^ (6) 

where the above equation is valid for all oscillators with 
odd A^. Performing the same type of analysis for even 
number of oscillators, where the nearest neighbors are 180 
degrees out of phase, we get: 

|z,|^-^"^" = (a + 3fc)'/' (7) 

Note that in the above equation the amplitude of the 
oscillation is independent of the total number of oscillators, 
while in Eqn. (6), this amplitude increases with increasing 

1/2 

N, asymptotically approaching (a -I- 3fc) as n — > oo. 

3. III. CONTRACTION ANALYSIS 

This section uses partial contraction analysis to analyze 
the stability of the synchronous and out-of-phase states 
for the system in Eqns. (1) - (3), for the A^ = 3 case. Here 
we use the partial contraction results first derived in the 
paper by Pham and Slotine, see Pham (2007). The results 
state a simple sufficient condition for global exponential 
stability on a flow-invariant linear subspace M (i.e. a linear 
subspace M such that Vt : F{M,t) C A^) as given by: 

-Klin (VLV"^) > SUpXmax ^ ) (8) 



where V forms a basis of the hnear subspace, A^^ (or- 
thogonal to Ai), dF/dx. is the Jacobian of the uncoupled 
system, with F given by Eqn. (2) and L is the coupling 
matrix, to be given below. The terms Xmin and Xmax 
indicate the minimum and maximum eigenvalues of the 
symmetric parts of the matrices L and dF/dx, respec- 
tively. Intuitively, the above condition insures that the 
system is contracting in the orthogonal subspace, A^"*", 
thereby insuring that the dynamics converge exponentially 
to M. 

Prom Eqn. (1), the coupling matrix, L for the N = 3 case 
is given by: 

/ / -/-/\ 
L = k -/ / -/ (9) 

V-/ -/ / / 

where I above is a 2 x 2 identity matrix. The entire phase 
space, A4 ® A4'^, is spanned by a total of six vectors: 
the two vectors spanning the synchronous solution, plus 
the four vectors spanning a linear vector subspace formed 
by the two degenerate out-of-phase solutions. For con- 
venience, let's define the subspace corresponding to the 
synchronous case as: 

Msync = {(x, X, x) : X e 7^^} (10) 

and the subspace corresponding to the out-of-phase states 
as: 

Mphase = {(^X, R^X, R^xj , (^X, R^ X, R^ X^ } (11) 

Where R is a 2 x 2 rotation matrix, with the rotation angles 
of {27r/3,47r/3} obtained from Eqn. (5) for the TV = 3 
case. The corresponding eigenvectors can be obtained by 
substituting x = {1, 0} and {0, 1} into Eqns. (10) and (11), 
resulting in six orthogonal eigenvectors. In complex form 
the two out-of-phase solutions given in Eqn. (11) can also 
be written as: z- (l, e'^''^^, e*'*'^/^) and z-{l, e'^'^^^, e^2,r/3^_ 
As before, the entire phase-space of solutions is spanned 
by the sum: Msync © Mphase- 

Having obtained all the eigenvectors, we can now use Eqn. 
(8) to analyze the stability of either the synchronous or the 
out-of-phase states, represented by Msync and Mphase, 
respectively. As will be shown shortly, this stability de- 
pends on the value of the coupling constant, k, with the 
synchronous solution being stable for negative k (or dif- 
fusive type of of coupling) and the out-of-phase subspace 
being stable for positive values of k (representing repulsive 
coupling) . 

To analyze the stability of the synchronous state, we 
equate: Msync = M and Mphase = M^, thereby ob- 
taining the 6x4 matrix V from Eqn. (11). The four 
eigenvalues of VLV'^, with L given in Eqn. (9), are all 
identical and given by: Ai.2,3.4 = 2/c. The eigenvalues of 
the symmetric part of the Jacobian, dF/dx, are given 
by: a — |xp and a — 3|a;p, which are upper-bounded by 
a. It follows from Eqn. (8) that the solution converges 
exponentially to Msync when, 

k < -a/2 (12) 

Next, doing the reverse by equating: Mphase = M and 
Msync = M-^, we again obtain the eigenvalues of VLV'^ 



(with V now given by Eqn. (10)): X^^ = —k. Combining 
the above results and again using Eqn. (8), we now get the 
condition for the global stability of Mphase- 

k>a (13) 

Equations (12) and (13) give conditions for the global 
stability of synchronous and out-of-phase dynamics, re- 
spectively. 

It is instructive to compare the model in Eqns. (1) and 
(2) to another coupling architecture which directly uses 
rotational matrices, see Seo (2007) and Pham (2007), to 
create globally stable out-of-phase state, in the following 
way (for a network of N oscillators): 

Xj=F(xj)+fc(xj-R2,Xj_i) (14) 



The above model results in a globally stable state where 
the nearest neighbors are displaced out-of-phase by 2tt/N. 
Unlike the model in Eqn. (14), the coupling term in Eqn. 
(1) does not need to be adjusted to get the 27r/A'' out- 
of-phase state as more oscillators are added, provided that 
the number, TV, is always increased by 2, so that N remains 
odd. As described in the previous section, for A^-even, the 
system splits into two identical synchronous groups 180 
degrees out-of-phase with each other. 

The other significant differences of the out-of-phase state 
created by the coupling term in Eqn. (14) from the system 
analyzed in the present work include: I. Existence of out- 
of-phase solution for any number, N, of oscillators. This 
is in contrast to the system analyzed here, where (as 
mentioned above), an odd N is needed for an out-of-phase 
state, II. The existence of a single (non-degenerate) out-of- 
phase solution, unlike the two solutions given in Eqn. (11), 

III. Different phase difference between nearest neighbors. 
Thus, the phase difference given by Eqn. (5) is such that 
the nearest neighbors are maximally out-of-phase, while 
the phase difference from rotational coupling is such that 
oscillator j is advanced from j — 1 by a phase of 2n/N, and 

IV. Preservation of uncoupled amplitude, ^/a in the out- 
of-phase state. This can be seen directly by substituting 
the out-of-phase solution xj = R2^Xj_i into Eqn. (14), 
whereby the coupling term drops out. 

Points I-III can be explained by pointing out that in Eqn. 
(14), the rotational coupling creates an out-of-phase state 
in the same way that a synchronous state is created in dif- 
fusive coupling. In other words, the oscillators all try to be 
in synch with the rotated by 2tt/N solution of their nearest 
neighbor. Therefore the out-of-phase state in rotational 
coupling is actually analogous to the synchronous state 
in diffusive coupling, which also has uniqueness, global 
stability (for appropriate values of k), and existence for 
any value of N. Point IV, that is the increase in the 
amplitude of oscillation (given by Eqn. (6)) in the out-of- 
phase state, is actually essential for the creation of Nui 
frequency oscillations when the two arrays are globally 
coupled. This model, inspired by the nervous system of 
a leech is analyzed in the following section. 



4. IV. GLOBALLY COUPLED ARRAYS AND THE 
ONSET OF HIGH FREQUENCY OSCILLATIONS 



Here we globally couple two arrays, each described by 
Eqns. (1) and (2), with the only difference being that the 
first array has repulsive coupling, given by fc,. > and the 
second array has diffusive coupling, given by kd < 0. The 
dynamics were inspired by the CPG-motorneuron network 
that controls the heartbeat of a medicinal leech, Buono 

(2004) . The heartbeat of the leech is driven by direct 
contact between two arrays of motorneurons, such that 
on one side of the leech the heart beats in a rear-to-front 
(peristaltic) fashion, well described by the out-of-phase 
state of coupled limit cycle oscillators. On the other side, 
the heart beats synchronously and is therefore represented 
by the diffusively coupled array, where the synchronous 
state is stable. The total system has the following form: 

N 

X- = F(xr) + kr (xj - xj+i - xj_i) -I- |x^| (15) 

k=l 

N 

xd = F(xf ) + k, (xf - x^Vi - xf_0 +cJ2 Ixlcl (16) 

k=l 

where as before, F(x) is defined in Eqn. (2). Based on the 
results of the previous section, we know that the system in 
Eqn. (15) has a stable out-of-phase state (corresponding 
to peristaltic motion), while the system in Eqn. (16) has 
a stable synchronous state. 

It has been shown both computationally (see Palacios 

(2005) ) and analytically (see Landsman (2006)) that a 
system of this type undergoes a bifurcation, as the global 
coupling constant, c, increases. For c > Cbif, the {x'*} 
array begins to oscillate in phase at the ultra-harmonic fre- 
quency given by Nlu. Following the method in Landsman 
and Schwartz, Landsman (2006), the bifurcation value of c 
that leads to high frequency oscillations can be calculated 
by solving for the value of the parameter P that causes a 
Hopf bifurcation in the following equation: 

Zj = {a ~ kd + iuj)zj — \zj\^Zj + P (17) 

where we have used the complex formulation, z — x + 
iy, of Eqn. (3) for the diffusively coupled array and 
substituted the synchronous solution: Zj = Zj+i = Zj-i. 
The bifurcation diagram for Eqn. (17) as a function of 
P is plotted in Figure (4), where bold lines at higher |P| 
indicate stable equilibria, with the broken line in the center 
showing an unstable equilibria. 

The bifurcation diagram in Fig. (4), in fact also corre- 
sponds to the bifurcation digram for the onset on high 
frequency oscillations in the diffusively coupled array, 
given by Eqn. (16). Namely the solid lines in the figure 
correspond to the Noj frequency oscillations in the x"^ 
array, while the broken line corresponds to oscillations 
close to the limit cycle frequency, uj in the same array. 
We can therefore obtain the bifurcation value of the global 
coupling, c for the onset of Nlo frequency oscillations by 
first solving for the bifurcation value of P (given by Pbif), 
and then using: P^if w CbifN\z^\ to calculate c^if. Here 
\z^\ is the amplitude of the out-of-phase oscillation in the 
repulsively coupled array. It is given by Eqn. (6), with 
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Fig. 3. NuUclines for Eqn. (17) for different values of P 
{a = 0.9, fed — —0.1, oj = 1/2). The green nuUclines 
are the i = nuUclines for various values of P. The 
red nuUcline is the y = nuUcline. Higher nuUclines 
corresponding to higher values of P. The circular orbit 
is a limit cycle for P = 0. 
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Fig. 4. Bifurcation diagram as a function of P, (a — 0.9, 
kd — —0.1, UJ = 1/2). Bold lines correspond to 
the high-frequency oscillations after the system in 
Eqn. (17) has been driven through a bifurcation for 
P>PMf. 

k ^ kr- Solving for Cf,i/, we obtain the value of the 
bifurcation to Nuj oscillations as a function of P: 

^ Pbif _ Pbif /-.„N 

""'f ^ N\z-^\ iV[a-f-fc^(H-2cos(7r/Ar))]i/2 ^ 

where Eqn. (6) was used in the denominator. 

The effect of the constant P on the dynamics of Eqn. 
(17) can be seen by referring to the nuUclines diagram, 
in Figure 3. The various green curves correspond to the 
i = nuUclines plotted for the corresponding values of 
P, with the P = nuUcline crossing the origin. Since 
P is a real constant in Eqn. (17), the y = nuUcline 



(shown in red) does not change with the parameter P. 
As |P| increases, the system bifurcates to a steady-state 
with a stable fixed point given by the intersection of the 
corresponding nullcHnes. The bifurcation value of \P\ is 
found at a point where the intersection of the vertical line 
with the X = nuUcline no longer has three real roots (see 
for example Guckenheimer (1983)), and given by; 
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where 7 = (a — kd) /3. Note that 7 > 0, since fc^ < in 
diffusive coupling. Substituting Eqn. (19) into Eqn. (18), 
we have the bifurcation value as a function of a, N, u, and 
the repulsive and diffusive coupling constants, given by kr 
and kd, respectively. 

The mechanism behind the onset of ultraharmonics for 
a similar type of coupling architecture was analyzed in 
Landsman (2006). Here we briefly summarize the mecha- 
nism behind the onset. The onset of high frequency oscil- 
lations hinges on the amplitude of the repulsively coupled 
array being higher than the amplitude of the diffusively 
coupled array. This results in a diffusively coupled array 
being driven through a bifurcation first, when c > c^if, 
with Ciiif given in Eqn. (18), and thereafter being driven by 
the out-of-phase dynamics of the repulsively coupled array, 
which causes the high-frequency synchronous oscillation 
in x'*. In the language of contraction theory, Lohmiller 
(1998), for c > ctif, the diffusively coupled array becomes 
a contracting system and can therefore be driven at the 
ultraharmonic frequency provided by the repulsively cou- 
pled array, which after the bifurcation acts like a drive. If 
c is increased even higher, beyond the bifurcation value of 
the repulsively coupled array, then total oscillator death 
in both arrays results. For c below this critical value, 
but above Ctif, given by Eqns. (18) (with P^if given by 
Eqn. (19)), high frequency synchronous oscillations are 
produced in the diffusively coupled array. 
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